Structure of Bi nanolines: using tight-binding to search parameter space 
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We describe how we have used tight binding calculations as a quick, efficient tool to search for 
possible structures of Bi nanolines on Si(001). After identifying promising candidate structures, we 
have concentrated on these with ab initio electronic structure techniques. The energetics of the 
tight binding are shown to be in good agreement with the density functional calculations and with 
experimental observations, and have proved invaluable in the search for a structure, validating the 
use of tightbinding as a search tool. 

PACS numbers: 



I. INTRODUCTION 

When confronted with a new structure observed in ex- 
periment, or when trying to find possible pathways for 
reactions or diffusion events, there is often a need for per- 
forming electronic structure calculations on many differ- 
ent structures rather quickly, in order to search through 
parameter space efficiently. However, the need is not 
just for speed, but also for reasonable accuracy, as pos- 
sible candidates must not be discarded out of hand. In 
this paper, we explore the use of tight binding as a tool 
to perform a search for the structure of Bi nanolines on 
Si(001), and consider future directions for these types of 
search using techniques developed over the last few years. 

If a Bi-covered Si(OOl) surface is annealed at 570- 
600° C, most of the Bi desorbs, but a small amount self- 
assembles torfarm nanolines, 1.5 nm wide and often over 
400 nm longufl. These .nanolines have great promise as 
templates for nanowiresO, when metals such as Al or Ga 
are deposited on them. However, to understand and con- 
trol these processes and the growth of the nanolines, it 
is important to know their, atomic structure; earlier pro- 
posals for the structureaErEl have been shown to be wrong 
by recent experimental results, and a valid structure has 
only recently been found using.a combination of STM and 
electronic structure modcllingtl We describe the search 
for the structure, and the tools used for that search, in 
this paper. n 

Tight bindingQ is a compromise between accuracy (it 
makes approximations, but preserves quantum mechan- 
ics) and computational efficiency (simulations of several 
hundred atoms can be run on a typical desktop work- 
station in a few hours). On semiconductor surfaces, 
where a rich variety of reconstructions exist, depend- 
ing on a subtle interplay of electronic and geometric ef- 
fects, the use of a quantum mechanical technique is vi- 
tal if correct answers are to be found. We have devel- 
oped parameterisations for Si(001), H on Si(001)@ and 
Bi on Si^OOjUO, and performed simulations cm-,lw|drpgeH| 
diffusioniiaO, gas-source growth of Si(001)B&Mliy 



and Bi on Si(001)U#Effl and have found that tight bind- 
ing produces results which are in good agreement with 
both experiment and first principles techniques. In this 
paper, we further explore the use of tight binding as a 
tool to search through parameter space. 

The problem we were faced with was that of finding 
the structure of the Bi nanolines. Various criteria had 
been found from experiment which possible candidates 
needed to fulfil: the lines must have a high kinking en- 
ergy; the lines must occupy four Si dimers but lie between 
the substrate dimers; they must be more stable than the 
Bi-(2xn) high coverage phase; and they must repel miss- 
ing dimer defects and steps down. We needed an accu- 
rate, quantum mechanical technique that would allow us 
to test many possible structures quickly, allow us to have 
confidence in new reconstructions and distorted bonds, 
and that would also allow us to model large unit cells 
(to test kinking and defect repulsion) - a linear scaling 
or O(N) tight binding technique was found to be per- 
fect. The rest of the paper describes the search: we start 
by describing the technical details of the modelling; we 
then consider the problem of modelling Bi on Si(001) and 
comparing structures with different amounts of Bi, and 
look at the stability of the Bi-(2xn) high coverage phase; 
we present the details of various more or less successful 
models, and compare their energies to first principles cal- 
culations, and then we conclude. 



II. TECHNICAL DETAILS 

We performed tight binding simulations using a linear 
scaling technique (where the computational effort scales 
linearly with the number of atoms., not with its cube as 
is true for traditional techniques^-!]) as large cells, of up 
to several thousand atoms, were necessary to test the 
kinking and dimer repulsion properties of the nanolinc 
structures. Specifically, we used apjmplcmcntaiiipn of 
the Density Matrix Method (DMM)E3 by GoringeN. We 
used tight binding parameterisations for Si-Si bonds de- 
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veloped specifically for Si(OOl)! and for Bi on Si(OOl)!. 
The unit cells were always ten layers deep, with the bot- 
tom two layers constrained to lie in bulk- like positions, 
and terminated in hydrogen. For simple structure cal- 
culations, each cell was one dimer row wide (using the 
p(2 x 2) reconstruction) and sixteen dimers long (i.e. 384 
atoms). When calculating kinking energies we used the 
same size cell, but four dimer rows wide (i.e. 1536 atoms). 
These methods have been used before with great success 
for calculations of many thousands of atomsj-Xe.g. for 
calculation of step kinking energies on Si(001)E3). 

For an accurate assessment of promising candidates, 
we performed first principles density functional (DFT) 
simulations within the locaJ density approximation 
(LDA) using the VASP codec3. We used ultrasoft pseu- 
dopotentials, a plane wave cutoff of 150 eV (which 
gave energy difference and force convergence) and a 
Monkhorst-PackEil grid of 4 x 4 x 1 k-points. Each unit 
cell was the same depth and width as the tight binding 
cells, but only eight dimers long, again with the bottom 
two layers held fixed in bulk-like positions and terminated 
in hydrogen. The change in energy caused by the shorter 
cell was found to be less than 0.1 eV; the longer tight 
binding cells were used to investigate nanoline-dimer re- 
pulsion, which could not be done using LDA. 

A key question to consider when modelling the stabil- 
ity of structures with different chemical elements is how 
to compare structures with different stoichiometries. As 
the experiments we are considering were performed at 
high temperature on a Si(001) substrate, we can use the 
energy of bulk silicon to correct for different numbers of 
silicon atoms (though this raises questions about k-point 
sampling and error cancellation for the LDA calculations 
that will be considered below). For bismuth, however, we 
are faced with a harder problem. In the next section, we 
detail the approach used for calculating the periodicity of 
the high coverage Bi-(2xn) phase, and extract a sensible 
reference energy for Bi dimers on Si(001). 



III. ENERGETICS OF BI ON SI(OOl) 

At high coverages, Bi forms rows of dimers on Si(001), 
relieving the strain due to the large difference in lattice 
constants by having missing dimer defects at regular in- 
tervals, forming a high coverage, Bi-(2xn) reconstruc- 
tion. As we see Bi nanolines forming within this (2xn) 
surface (and we also see that the nanolines are more sta- 
ble than it - they form while it is desorbing) we start 
by describing it. We approach this problem, which will 
involve different amounts of Bi, and hence different sto- 
ichiometrics, in the same way that wa-have approached 
modelling of the Gc/Si(001) surfaceEl It is perfectly 
easy to calculate the energy of formation of trenches of 
missing dimers simply by removing Bi dimers, but to 
make contact with experimental results we need to find 
an energy for these Bi dimers that we have removed. We 
shall assume that the Bi-covered surface is in thermal 



equilibrium (which is easy to arrange experimentally), so 
that for a given number of Bi atoms we can write the 
probability of any given arrangement of these atoms as 
proportional to exp(E a /kT), where E a is the energy of 
the arrangement. The assumption of thermal equilib- 
rium greatly simplifies the system that we are consid- 
ering, while being perfectly legitimate at a number of 
different temperatures. 

Now, as we are in thermal equilibrium, we only need 
know the relative energetics of different arrangements of 
the given number of Bi atoms. It will be useful to de- 
fine a reference arrangement which we take to be a non- 
defective, lattice-matched layer of Bi covering a certain 
area of the Si(001) surface (with a boundary which we 
assume to be rectangular, but whose shape is arbitrary). 
Now we consider arrangements where we form regularly 
spaced trenches of missing dimers; conceptually we start 
by adding a number of dimers to the reference system 
equal to the number that we will remove and then form- 
ing the missing dimer trenches and relaxing. This con- 
ceptual approach allows us to compare different spacings 
by using the energy of Bi dimers in the perfect mono- 
layer to correct for different amounts of Bi, which we 
find by calculating the difference in energy between a 
clean Si(001) surface and the same surface covered by a 
monolayer of Bi. Using our tight binding parameterisa- 
tion we find this number to be -11.41 eV per Bi dimer. 

We show the difference in energy per Bi dimer between 
a full, strained monolayer and the reconstructed surface 
for various values of period, n, in Fig. |. There is a well 
defined minimum for periods of four and five (i.e. three 
and four Bi dimers) whichJsJn excellent agreement with 
experimental observationarEj. The modelling of this re- 
construction gives us further confidence in our parame- 
terisation, and presents a limit for the energetics of the Bi 
nanolinc (which must be more stable than the Bi-(2xn) 
surface). Our tight binding simulations find that the en- 
ergy per Bi dimer (adsorption energy and surface energy) 
in this reconstruction is -11.6 eV at n= 4 and n= 5. 

When considering structures with different amounts of 
Bi, we can use the fact that the nanolines are observed 
to form in the high coverage Bi-(2xn) phase to extract 
a sensible reference energy for Bi dimers. We have cho- 
sen the energy to be -11.4 eV per dimer, the energy of 
the reference phase, and used this where necessary. As 
mentioned above, we use the bulk Si energy to correct 
for differing numbers of Si atoms. 



IV. DETAILS OF SUCCESSFUL MODELS 

In this section, we present details of various of the 
models considered, and discuss how well tight binding 
has performed as a tool to discriminate between these 
models in the light of subsequent DFT modelling. First, 
we start by summarising the requirements placed on the 
models by experimental observations. 




Periodicity (dimers) 



FIG. 1: Energy difference per dimer between reconstructed 
surface and full monolayer for Bi on Si(OOl) in eV, plot- 
ted against periodicity of reconstruction (where number of 
Bi dimers is one less than periodicity). 



A. Experimental Constraints 

We have been able to distil a number of requumuents 
for the nanoline structure from our experimentsBEItm: 

• Stability: The lines form at a temperature where 
the high coverage (2xn) Bi islands are evaporating, 
which implies that the energy per Bi dimer must 
be larger than these islands (in tight binding, more 
than -11.6 eV/dimer as given above). In RHEED 
experiments, the difference in desorption barrier for 
Bi from the islands ancLfrom the nanolines was de- 
termined to be 0.25 eVtl. 

• Width, Registration and Appearance: The nano- 
lines run perpendicular to the Si(001) surface 
dimer rows, occupying the space of four Si dimers 
(15.36 A). They consist of two dimer-like features 
separated by ~6.3 Awhich lie between the Si(001) 
dimers, but are in line with the dimer rows. 

• Straightness: Over many hundred observations of 
many lines (for instance in previous paperso'B), a 
kink has never been seen. Given the high tem- 
perature at which these nanolines form, it seems 
likely that a large thermodynamic kink energy is 
required, though a kinetic barrier cannot be ruled 
out. 

• Defect Exclusion Zone: The nanolines are 
observedO to repel defects and down step edges to 
a distance of 3-4 nm. As the defects and steps have 
a tensile stress field, we expect that the nanoline 
model should also have a tensile stress field. 
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FIG. 2: Ball-and-stick models of various candidate modaJe 
for the Binnanoline. Labels, from top to bottom: originalcl, 
TwolDVsEI, B-type 4DV, two ad-dimers and two squares. En- 
ergetics are given in table |, and the structures are discussed 
more fully in the text. 



B. Energetics from tight binding and DFT 

Over the course of approximately a year, we have 
considered many different models. In general, we have 
first calculated the energy and kinking energy of possible 
models using tight binding. Only if these looked promis- 
ing were they modelled with the more expensive DFT. 

Within this section we present details of five unsuccess- 
ful models and what we believe to be the correct modeltj, 
giving energies found with both tight binding and DFT. 
These results show that the tight binding calculations 
certainly ordered the structures correctly in terms of en- 
ergy compared to DFT calculations (to within 0.1 eV) 
even though the absolute energies were sometimes rather 
too large, and thus confirm the use of this method as 
extremely valuable when searching a large configuration 
space. It is also worth noting that the tight binding simu- 
lations were computationally cheap (running in a matter 
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FIG. 3: Ball-and-stick model for the successful modelQ. There 
has been significant reconstruction below the Bi ad-dimers, 
yielding 5- and 7-membered rings and changes down to the 
fifth layer. However, this structure is extremely stable in both 
TB and DFT, and satisfies all experimental constraints. 



of hours on a desktop workstation) and also that they 
deepened the collaboration between experiment and the- 
ory (many tight binding simulations were run by an ex- 
perimentalist (JHGO), enabling detailed discussions of 
the remaining promising candidates). As an illustration, 
we found that tight binding relaxation calculations on a 
large unit cell used for kinking calculation (1536 atoms) 
took 2 hours and 4 minutes on a PowerMac G4/867MHz 
(using open-source g77 and ATLAS BLAS routines op- 
timised for AltiVec) and 1 hour and 45 minutes on a 
Pentium 3/800MHz (using the Intel Fortran compiler 
for Linux and either ATLAS BLAS routines or the In- 
tel Math Kernel library). 

In Fig. H we show the ball-and-stick models of five un- 
successful-candidates, and in Fig. |3| we show the final 
candidates (shown in off-axis views as the structures are 
rather complex) . The five unsuccessful models shown are: 
the original structure proposedBBo, which fitted original 
STM data but is only 3 Si dimers wide and has a very low 
kinking energy (labelled "Original" ) ; the model proposed 
by Naitoh et alu, which has the wrong spacing between 
features, is energetically poor and has a very low kink- 
ing energy (labelled "TwolDVs"); a structure made from 
back-to-back rebonded B-type steps, with Bi substituted 
in the second layer, which has the wrong appearance and 
is not energetically good enough (labelled "Btype4DV" ) ; 
two Bi ad-dimers, which have no kinking energy and are 
energetically poor, and have no incentive to remain in 
this configuration (labelled "TwoAdDimers"); and two 
squares of Bi-pn top of a reconstruction similar to the fi- 
nal structureu and to the structure of steps on As-covered 
Ge(001)EI 

We give the energies for these models, along with a 
variation on the Btype4DV structure (the same struc- 
ture but with only a 3 dimer width, which is energet- 
ically rather better, but has the wrong spacing and a 
low kinking energy, labelled "Btype3DV" and not illus- 
trated) and the final candidate (labelled "DoubleCore" ) 
in Table |. Where the number of silicon atoms has dif- 
fered from the clean surface, we have corrected using an 
appropriate value for the bulk energy for silicon. In the 
DFT calculations, it is questionable to use simply an en- 



TABLE I: Energies of various candidate structures as found 
by tight binding and DFT. Energies arengiven relative to the 
three dimer model originally proposedBB, but now known to 
be inconsistent with experimental observations. Each model 
contains two Bi dimers, and has been corrected for differ- 
ent numbers of Si atoms as described in the text, except for 
TwoSquares, where the energy is given per square (i.e. equiv- 
alent to the others). Labels are explained in the text, and 
models are shown in Figs. U and pi 



Model TB (eV) DFT (eV) 

TwoSquares 1.179 1.363 

TwoAdDimers 1.565 0.848 

TwolDVs 1.232 0.759 

Btype4DV 0.611 0.564 

Btype3DV 0.473 0.279 

Original 0.000 0.000 

DoubleCore -1.904 -0.720 



ergy calculated from a bulk Si cell, as the k-point sam- 
pling will differ unless fully converged in both cells; we 
also expect error cancellation to be better in the same 
sized unit cells. We approached the problem from two 
directions: first, we calculated the total energy of a clean 
Si(001) surface using both an eight layer and a ten layer 
slab, and used the energy difference to give us a value for 
the energy of bulk Si. We compared this to calculations 
of a bulk Si cell which were fully converged with respect 
to k-point sampling (we used a 4 x 4 x 4 mesh, and 
checked with an 8 x 8 x 8 mesh), and found that the 
energy difference less than 0.05 eV per atom. The energy 
that we used was 5.96 eV/atom. 

From these energies, we can see that the tight bind- 
ing calculations have correctly ordered almost all of the 
structures, though the energies are sometimes too large. 
It is extremely pleasing that it has performed so well 
for such a large variety of structures, especially the final 
candidate which contains a large stress and significant re- 
construction. In particular, it can be seen that there are 
a wide range of bond lengths (Bi-Bi are generally about 
2.9 A, but Bi-Si vary from 2.60 to 2.82 A) and angles 
(from 84° to 106°) which the tight binding calculations 
have modelled correctly. Bi seems to favour p 3 bonding, 
with 90° bond angles and a pair of electrons in an s or- 
bital - and the tight binding calculations have modelled 
this well. The only point where the ordering is wrong is 
for the TwoSquares structure, which tight binding found 
to be more stable than it should have been. The Bi-Bi 
bonds are generally shorter than found in DFT calcula- 
tions (~3. 0-3.1 A), suggesting that the Bi-Bi bond in TB 
might be too strong, accounting for the increased stabil- 
ity of this structure. The magnitudes of the TB energies 
are too large, probably because the model was fitted to 
bulk Bi-Bi and Bi-Si bonds, and this environment is at 
the surface, allowing Bi to form bonds with angles close 
to 90° (its preferred bonding angle). 

The unsuccessful models were rejected for a num- 
ber of reasons: they were energetically poor (TwolDVs, 
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TwoAdDimers, TwoSquares); they had the wrong ap- 
pearance (Btype4DV and Btype3DV); they had a low 
kinking energy (Original, TwolDVs, TwoAdDimers, 
Btype3DV); they had the wrong spacing of Bi dimers 
(TwolDVs, Btype3DV, Btype4DV); they had the wrong 
registry of Bi dimers relative to the substrate (Origi- 
nal, TwolDVs); and they were less stable than the high 
coverage Bi-(2xn) surface (TwolDVs, TwoAdDimers, 
Btype3DV, Btype4DV, TwoSquares). 

The final candidate (DoubleCore) has the correct ap- 
pearance, a high kinking energy (3.75 eV from tight bind- 
ing calculations), the correct spacing of features (6.3 A), 
a high repulsion for defects and down step edges (out to 
6 or 8 dimers or ^3 nm), and is more stable than the 
high coverage Bi-(2xn) surface. We have discussed this 
structure in more detail elsewhereH. 



V. DISCUSSION AND CONCLUSIONS 

The key drawback with tight binding is that a parame- 
terisation must be created for every interaction, and this 
is frequently a time-consuming procedure. Another ex- 
tremely important problem for tight binding is that of 
transferrability: how well a parameterisation describes 
interactions that are dissimilar to the environment used 
for fitting. The Bi-Bi and Bi^3i parameterisation used 
was created simply and quicklyO, using the universaLtight 
binding parameters of Harrison for Bi-Bi hoppingcj and 
fitting the scaling to the bulk modulus only; the parame- 
terisations have only been used in environments close to 
those used for fitting. This simple prpcedure (which has 
also been applied to Ge-Ge and Ge-SiH3) has been shown 
here to be remarkably effective. 

However, there is clearly a need for more reliable ways 
to explore parameter space in situations such as that en- 
countered in this paper, and developments in the field 
known as ab initio tight binding are providin g th e an- 
swers. Starting with the Harris-Foulkes theoryc3E3, the 
link between DFT methods using local, atomic-like or- 
bitals and tight binding was shown. This was quan- 
tified in an elegant demonstration that tight binding 
could be derived, from DFT via a well-defined series of 
approximations^. Ideas such as these have led to differ- 
ent approaches to improving and generalising tight bind- 
ing. 

Typically, confined atomic orbitals (found by solving 
the Schrodinger equation within a confining potential)E2l 



are used as a basis set, adding orbitals beyond the mini- 
mal set as necessary for completeness. The problem with 
any approach like this is to decide where to make ap- 
proximations. There are two key strands to these ap- 
proximations: to neglect self consistency (but to retain 
three-centre integrals)E3; and Xa add self-consistency but 
to neglect three-centre termsEi Of course, there have 
been developments such as adding self-consistency to the 
three-centre approached 

The ultimate generalisation of these ideas lies j-ki 
local oxhital-based DFT methods such as SiestaES, 
PlatoEJH and Conquesto, which point the way for- 
ward, as they enable the creation of the matrix elements 
(equivalent to a parameterisation) accurately, both in 
advance and in real time: the link between tight bind- 
ing methods and DFT-based methods is becoming ever 
stronger. Clearly the level of approximation chosen needs 
to trade off against the size and complexity of the system 
being modelled. 

In conclusion, we have presented tight binding and 
DFT calculations for the high coverage Bi-(2xn) sur- 
face which are in good agreement with experiment, and 
we have explored a variety of possible structures for Bi 
nanolines. Tight binding calculations have been shown 
to agree well with DFT calculations in terms of order- 
ing structures, and enabled us to perform quick, efficient 
searches of parameter space. In particular, it has enabled 
us to find an extraordinary structure for the Bi nanoline, 
which in turn sheds light on the role of surface stress 
in self-assembly of nanostructures on semiconductor sur- 
faces. 
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